Longitudinal variation in resilient psychosocial functioning is associated with ongoing cortical myelination and functional reorganization during adolescence

Adolescence is a period of dynamic brain remodeling and susceptibility to psychiatric risk factors, mediated by the protracted consolidation of association cortices. Here, we investigated whether longitudinal variation in adolescents’ resilience to psychosocial stressors during this vulnerable period is associated with ongoing myeloarchitectural maturation and consolidation of functional networks. We used repeated myelin-sensitive Magnetic Transfer (MT) and resting-state functional neuroimaging (n = 141), and captured adversity exposure by adverse life events, dysfunctional family settings, and socio-economic status at two timepoints, one to two years apart. Development toward more resilient psychosocial functioning was associated with increasing myelination in the anterolateral prefrontal cortex, which showed stabilized functional connectivity. Studying depth-specific intracortical MT profiles and the cortex-wide synchronization of myeloarchitectural maturation, we further observed wide-spread myeloarchitectural reconfiguration of association cortices paralleled by attenuated functional reorganization with increasingly resilient outcomes. Together, resilient/susceptible psychosocial functioning showed considerable intra-individual change associated with multi-modal cortical refinement processes at the local and system-level.


Distribution of distress scores and resilient psychosocial functioning with respect to demographics
Having derived distress scores (see Supplementary Table S1 for loadings) and resilient psychosocial functioning (ResPSF) scores and their intra-individual rates of change (∆), we tested whether scores differed between self-reported sexes or as an effect of age (see Supplementary Table S2 for statistics).In the larger sample of 712 individuals, we observed significantly higher ResPSF and lower distress levels in males than in females.This sex difference was not significant in the smaller sample (n = 141) included in the imaging analyses.Neither ResPSF nor distress scores were significantly associated with age.Intraindividual change in ResPSF and distress showed no sex or age effects in neither sub-sample.As change in ResPSF was the main behavioral variable of interest in the current study and did not show a significant sex difference, we did not perform further sex-stratified analyses.

Potential links between changes in ResPSF and changes in adversity exposure
Given that ResPSF were computed separately for each time point by predicting distress levels based on adversity levels measured at each time point, our approach inherently adjusts for changes in adversity exposure between measurement time points.We still wanted to confirm that changes in ResPSF were not simply a correlate of increasing or decreasing adversity levels.Correlating ∆ResPSF with ∆s of individual risk assessment scores revealed no significant associations between ∆ResPSF and ∆-scores from the MOPS (r = 0.10), APQ (r = 0.10), LEQ (r = -0.02;all p> 0.05), or SES (remained unchanged for most individuals; Supplementary Figure S1).Changes in adversity exposure were further not associated with changes in myelin-sensitive MT (all p> 0.05).

Supplementary Figure S1. Changes in adversity exposure between measurement timepoints. The upper row depicts intra-individual changes in adversity measures (n = 141 individuals). The lower row depicts respective
Pearson's correlations between changes in adversity measures and changes in resilient psychosocial functioning (∆ResPSF; n=141).No significant association with ∆ResPSF was observed for change in adversity exposure as captured by the Measure of Parenting Style (MOPS; r = 0.10, p = 0.24), the Alabama parenting questionnaire (APQ; r = 0.10, p = 0.21), or the Life events questionnaire (LEQ; r = -0.02,p = 0.79).IMD = Index of mean deprivation.

Representativeness of the MRI subsample
To address concerns that developmental neuroimaging studies may not be representative of the general developmental population 1 , we descriptively compared the distributions of levels of environmental risk exposures as well as behavioral outcome measures in n = 144 individuals included in neuroimaging analyses (the imaging subsample) and the larger sample of individuals included in behavioral analyses only (the non-imaging subsample; Supplementary Figure S2).We generally observed strongly overlapping distributions for distress and resilient psychosocial functioning (ResPSF) scores, as well as for APQ, MOPS, LEQ and CTQ questionnaires, between the two sub-samples.For SES, considered here as an index of mean deprivation (IMD), we observed a comparable range but a relatively higher proportion of higher SES in the imaging sub-sample (i.e., a more left-skewed distribution).This suggests a potential oversampling of individuals from higher socioeconomic backgrounds for the neuroimaging analysis.At the same time, it should be noted that the behavioral outcome measure (ResPSF) was computed in the full sample, including more individuals with lower SES, and also that SES was weighted with the lowest feature importance by the model used to predict ResPSF (Figure 1A).
We also note that the NSPN sample is a locally collected sample from London and Cambridgeshire in the UK.75% of the total sample and 84% of the imaging sub-sample were white.Moreover, the sample is predominantly healthy, thereby excluding individuals with, for example, neurodevelopmental disorders who are part of the general population.It will be important to test the replicability of the current results in more socio-economically and ethnically representative samples of the general population.The efficacy of potential resilience factors implied by the current study in a largely white, relatively affluent and healthy sample requires further validation in ethnically, socioeconomically and clinically defined groups that are underrepresented in this sample.

Sensitivity tests for ∆ResPSF * ∆MT effects
We tested whether the observed association between ∆ResPSF and ∆MT was robust to analytical choices by spatially correlating (Pearson's r) the original unthresholded t-map with t-maps derived from alternative analytical approaches (Supplementary Figure S3A).We observed virtually unchanged results when 1) not including mean ResPSF as a covariate in the general model (r = 1 with the original map), 2) including a quadratic age term to control for non-linear age effects (r = 0.998 with the original map), and 3) not winsorizing the input data (correlation of r = 0.987 with the original map).Next, we tested whether the observed prefrontal effect was present in different sub-samples.For this purpose, we drew 1000 sub-samples each containing 80% of the individuals and repeated the analysis 1000 times.The average correlation between the original t-map and the t-maps derived from the sub-samples was r = 0.929 (Supplementary Figure S3B).

Supplementary Figure S3. Sensitivity tests for ∆ResPSF * ∆MT effects. A)
Pearson's correlation between original the t-map (x-axis) presented in the main manuscript and t-maps derived from alternative analytical approaches (n = 141 individuals, all p < 0.001).B) Histogram depicting Pearson's correlations between the original t-map and t-maps computed based on 1000 randomly drawn sub-samples (80% of data).

Organizational axis of developmental change in MT
To assess whether ∆MT occurred in a synchronized manner across the cortex, we further observed that inter-regional covariance of ∆MT was organized along an anterior to posterior pattern (Figure S4).We then examined whether the topological distribution of ∆ ResPSF effects on ∆MT aligns with general organizational principles of intra-individual MT development and observed a positive correlation between the t-map and the previously identified principal axis of ∆MT (r = 0.46, pspin = 0.016, CI = [0.36,0.53]).Significance was assessed by a spin test (10000 permutations) correcting for spatial autocorrelations between cortical maps 2 .Thus, the association between change in resilient psychosocial functioning and MT maturation follows general organizational principles of MT development.This suggests that differences in the cortex-wide embedding of developmental change in anterior vs. posterior regions are reflected in the degree to which both apexes play a role in the development of resilient psychosocial functioning.The observed axis likely reflects differences in the timeline on which regions are most developmentally active with respect to myelination, and thus have differential relevance for adolescent resilience.

Cross-sectional effects of ResPSF and ∆ ResPSF
We observed no cross-sectional association between baseline MT and baseline ResPSF.However, we observed a negative effect of ∆ ResPSF on baseline MT in a parcel in the left middle frontal gyrus (L_p10p; t = -4.16,p10,000 permutations + FDR < 0.05).Thus, individuals who had lower ResPSF at baseline compared to the follow-up time point also had cross-sectionally lower levels of MT in this region (Figure S5).

PFC sub-clusters in ∆FC analyses
Using functional connectivity data, we investigated whether regions that show differences in ∆MT as a function of ∆ResPSF also exhibit different functional embedding.For this analysis, presented in the main manuscript, connectivity profiles were averaged across parcels that are part of the prefrontal ROI within each participant.As the mask/ROI includes large parts of the PFC and thus a nexus of different functional profiles, we further subdivided the cluster based on parcels' assignment to functional communities 3 .The ROI spanned sub-regions of frontoparietal, default mode, and limbic networks.Because large parts of the limbic network were excluded in FC analyses in this study due to low signalto-noise ratios (see Methods), we tested associations between ∆ResPSF and ∆FC only for frontoparietal and default mode sub-regions.
The prefrontal cluster showed predominantly negative change in FC across large parts of the frontal, temporal, and occipital cortex, as well as the insula and anterior cingulate cortex (Supplementary Figure S6A).Studying the association between ∆ResPSF and ∆FC in sub-parts of the prefrontal cluster, it appeared that positive associations were largely driven by sub-regions that are part of the default mode network.That is, when defining sub-regions that are part of the frontoparietal network as seed, two regions (right medial PFC and PCC; Supplementary Figure S6B) showed significant effects after 10.000 permutations and FDR (alpha = 0.05).In contrast, defining sub-regions that are part of the default mode network as seed revealed positive associations between ∆ResPSF and ∆FC in 21 regions, spanning medial and lateral PFC, posterior cingulate cortex, and parts of unimodal sensorimotor cortices (Supplementary Figure S6C).

Supplementary Figure S6. Effects of change in resilient psychosocial functioning (∆ResPSF) scores on change in functional connectivity (FC) of sub-parts of the prefrontal region-of interest presented in the main analysis (n=141 individuals). A) Group-average pattern of change in FC (∆FC) for the total prefrontal region of interest (ROI). B)
We observed no significant association between ∆ResPSF and ∆FC in the ROI sub-regions that are part of the frontoparietal network (FPN).C) Associations between ∆ResPSF and ∆FC in ROI sub-regions that are part of the default mode network (DMN).B and C include age, sex, site, and mean ResPSF as covariates.The respective ROIs are masked in black.Regions excluded due to low signal-to-noise ratio are masked in dark grey.

Robustness, estimated maturational processes, and future utility of the MPC maturational index
The maturational index based on microstructural profile similarity (MIMPC ) is a structural extension of the previously established functional connectivity maturational index, capturing conservative (i.e.strengthening of existing connections) and disruptive maturational modes (i.e., a reorganization of existing connectivity patterns 4 ).While the functional maturational index identifies increasing levels of reorganization from unimodal (little reorganization) to transmodal cortex (more reorganization), the microstructural maturational index mirrors insights gained from a previously established cortical topology of synchronized age effects on microstructural profile covariance 5 .In the MIMPC, we observe the strongest re-organization in frontoparietal heteromodal cortex, whereas a "frame" of ventral/paralimbic and dorsal/somatosensory cortex shows mostly an age-related strengthening of existing MPC patterns (i.e., little re-organization).This is consistent with what Paquola & Bethlehem et al. reported: association cortical areas, in which overall intra-cortical myelin content increases, develop towards a more "sensory" architecture, whereas regions in which preferably mid-to-deeper layers show increases in myelin develop towards a more "paralimbic" architecture.This differentiation process is reflected in the 'disruptive re-organization' captured by the MIMPC, and mirrors the modular segregation observed in tractography-based adolescent data 6 .At the same time, paralimbic-temporal/ventral and somatosensory/dorsal regions show 'conservative development' (i.e., little re-organization) in the MIMPC, suggesting that MPC patterns are well-defined prior to adolescence 5,7 .Overall, the MIMPC pattern meaningfully captures synchronized microstructural maturation and re-organization, consistent with previous observations.
To probe the robustness of the MI derived from microstructural data, we repeated the computation of the MIMPC based on different subsamples.First, we drew 100 sub-samples, each containing 80% of the individuals per NSPN age bin and repeated the analysis 100 times.The average correlation between the MIMPC map based on all individuals for whom MT data was available (n = 295 subjects, 512 sessions/datapoints) and MIMPC maps derived from 80%-sub-samples was r = 0.96 (Supplementary Figure S7A).Next, we assessed whether the MIMPC map can be observed with smaller sample sizes.To this end, we again drew sub-samples per age bin, but this time the size of the subsamples ranged between 20% and 100% (in steps of 5%) of individuals (see Supplementary Figure S7B).Last, results stayed consistent when computing the MIMPC -which reflects the correlation between baseline and age-related change patterns -based on Spearman's or Pearson's correlation (r = 0.99; pspin < .0001)Supplementary Figure S7C).
The maturational index based on microstructural profile covariance adds to our understanding of inter-regionally synchronized cortical maturation, capturing adolescent re-organization (integration and segregation) of primarily frontoparietal association cortex.That is, emerging work highlights integrated multi-scale approaches to elucidate biological risk factors associated with neuropsychiatric conditions.It is increasingly recognized that pathological functional perturbations are coupled with microstructural perturbations 8-11 .Moreover, taking a nuanced approach to studying intracortical myeloarchitectural profiles, beyond mean myelin content, has revealed parallel maturational processes at different scales and topologies 5,12,13 .The use of the MIMPC in future studies may mirror the current use of similar, already established measures such as the main axis of MPC age effects and the maturational index for functional connectivity.That is, the main axis of MPC age effects has already been linked to the cortical topology from microstructural profiles and histology 5 , which in tun has been combined with other measures of 'cortical wiring' in adolescence 14 .Assessing the topology of synchronized structural maturation based on intra-cortical profiles further extends classical structural covariance approaches, assessing e.g.cortical thickness covariance 15 , towards a more nuanced myeloarchitecture.Last, similar to our study, the previously established maturational index of functional connectivity 4 has been demonstrated to robustly capture sex differences in adolescent functional network maturation 16 .In summary, the MIMPC may be of interest for future studies of adolescent cortical maturation.

Supplementary Figure S7. Robustness of the MPC maturational index. A)
Histogram depicting correlations between the microstructural profile covariance maturational index (MIMPC) derived from all n = 295 individuals (512 sessions) and the MIMPC based on 100 sub-samples (80% of data).B) Correlations between the MIMPC derived from all 295 individuals (512 sessions) and the MIMPC based on subsamples of different sizes, ranging from 25-100% in steps of 5%.C) The MI is computed as the correlation between baseline and change patterns.Here, we applied both Spearman's and Pearson's correlation to compute the MI and correlated the MI pattern resulting from either method via a Pearson's correlation (r = 0.99, pspin <.0001).

Conservative vs. disruptive trends MIMPC Group differences
We tested for differences in the maturational index (MI) between groups of individuals who showed increasingly resilient vs. susceptible outcomes.Table S3 summarizes all 43 parcels displaying significant group differences in MIMPC (as determined both via z-tests and non-parametric permutation testing).Trends were determined with respect to the full sample MIMPC presented in Figure 2.For example, a positive ∆MIMPC in a region labeled 'conservative' in the full sample MIMPC is interpreted as '+∆ResPSF more conservative', whereas a positive ∆MIMPC in a region labeled 'disruptive' in the full sample is interpreted as '+ResPSF less disruptive'.Regions that show a significant group difference but neither a clear conservative nor disruptive pattern in the full sample are labeled as 'tipping points'.
Table S3.ROIs with significant group differences in MIMPC and their trends.Group differences were derived from z-tests and adjusted for multiple comparisons by thresholding at pFDR <0.05 as well as non-parametric permutation testing using 10.000 permutations.Tests were two-sided.

Supplementary Discussion
Insights gained and potential sources of noise related to using residuals as resilience scores Assessing resilient/susceptible outcomes based on the deviation of observed from expected levels of well-being given a certain adversity exposure is a well-established approach in the resilience literature [18][19][20][21][22][23] .It addresses the issue that simple quantification of mental health variables can provide only limited insights into resilience or susceptibility, as they are heavily conflated with individual differences in adversity exposure (Kalisch et al., 2021).Leveraging residuals to quantify better or worse well-being than predicted by adversity exposure thus provides a corrected well-being score that is adjusted for individual differences in stressor exposure and thus allows comparing resilience scores of individuals with differing exposure levels.That is, 'residuals' in the resilience use case are interpreted as 'residual variance in mental health problems' that is not explained by the normative response to exposure, and therefore indicate individually weaker response (resilience) or a stronger response (susceptibility).It is thus analogous to the process of correcting a dependent variable for potential confounders such as age or sex.
While the basic assumption of this approach is that residualized mental health outcomes reflect degrees of resilience/susceptibility, there are other potential influences such as 1) measurement error and noise related to the questionnaires used, 2) noise related to the performance of our prediction model, and 3) confounding influences of other environmental influences.
1) Self-report / retrospectivity biases may pose one source of measurement error, e.g. for the reporting of childhood maltreatment 25 , despite the relatively short reporting time windows of the current study.It should be noted, however, that such sources of measurement error are not unique to / caused by the residual approach, but rather persist from the original measures of psychosocial well-being.The resilience scores should therefore not contain more measurement error than the original measures of psychosocial well-being.In the current study, our longitudinal approach may further mitigate some aspects of measurement error if it is linked to e.g., self-report bias.That is, an individual systematically self-reporting his/her well-being as a bit better than it is will receive higher resilience scores at all time points.However, as our study investigates intra-individual change, such a general offset in resilience scores would not affect change scores, whereas it would affect crosssectional analyses more strongly.2) Another question is how well our prediction model controls for differences in adversity exposure in order to reveal resilient / susceptible responses -and to what extend influence is still uncontrolled for, thus causing noise in in our measure of individual resilience.The amount of variance explained by our model (approx.21%) is comparable with previous work reporting 24% of variance explained for the prediction of psychosocial functioning from family experiences in longitudinal settings 23 ; 21% when predicting internalizing symptoms from general life stressors during the covid pandemic 26 ; or 28% when predicting psychosocial functioning from childhood adversity 27 .We thus conclude that our model controls for exposure to a comparable degree as common resilience models.3) Capturing meaningful variation in psychological outcomes is complicated by the complex influence of a multitude of interacting factors, which likely contribute to the variance not explained by our model.Such factors may include genetic predispositions, other environmental risk or protective factors not measured here, but likely also sources of noise we cannot quantify.For instance, an individual with a genetic predisposition for mental illness may systematically show lower psychological well-being, which cannot fully be explained by adversity exposure and would thus create a bias in that individuals' derived resilience scores.Similar to 1), if general offsets in resilience scores exist due to e.g., a genetic predisposition, studying longitudinal change helps us to account for such an offset.
Overall, our results using resilient psychosocial functioning scores suggest a central role for multi-modal prefrontal maturation and more wide-spread re-organization of association cortices for resilience and susceptibility during adolescence, tested against null models by non-parametric permutation.This observation is well in line with previous reports of structural and functional involvement of these brain regions in stress responses and susceptibility/resilience [28][29][30][31][32] .We believe the fact that observed associations of residuals / resilience scores with multi-modal measures of cortical maturation survived non-parametric permutation tests, and were found in regions previously suggested by the literature, argues for a dominance of meaningful variance reflected in the scores used.
The Warwick-Edinburgh Mental Wellbeing Scale (WEMWBS; 39 ) is a 14-item questionnaire that measures mental well-being in a positively worded fashion.Participants were asked to respond on a 5-point Likert scale ('none of the time', 'rarely', 'some of the time', 'often', 'all of the time') to which degree statements described their experiences in the last two weeks.Example items are 'I've been feeling good about myself', 'I've been feeling useful', I've been feeling close to other people'.A sum score was used for present analyses.
The Schizotypal Personality Questionnaire (SPQ; 40 ) originally includes 74 binary ('present' or 'absent') self-report items designed to capture symptoms associated with the DSM-III definition of Schizotypal Personality disorder, such as psychotic-like experiences.We included only items that have previously been tested to be significantly associated and showing medium to high effect sizes with psychotic-like experiences as measured by the semi-structured PLIKS interview (PLIKSi; total score, hallucinations, delusions, and perceptual abnormalities; 35,41 .Items retained based on this face validity include: SPQ 4, 9, 13, 28,31, 40, 55, 60, 61, 63 and 64.

Adversity questionnaires information
The Life Events Questionnaire (LEQ; 42 asks individuals about significant life events that have occurred within the previous 18 months.Such significant life events include changing schools / college / jobs, moving, changes in family composition like death or divorce, disasters at home (e.g.fire), serious illness and/or hospitalization of self or someone in the close network of family and friends, deaths, loss of family pet, problems with or end of friendships, and others.Participants were instructed to focus on the most impactful event if they experienced multiple situations captured by the same category.Moreover, participants rated how (un)pleasant the event was ('very pleasant', 'pleasant', 'neither', 'quite unpleasant', 'very unpleasant') and whether it impacted them for more than 2 weeks.For our analyses, we used the LEQ sum score capturing how many adverse life events (i.e., events scored as 'quite unpleasant' or 'very unpleasant') an individual faced within the given time period.
The Child Trauma Questionnaire (CTQ; 43 ) measures abuse and neglect up to age 18 years on a five-point Likert scale ('never' to 'always') in five overarching categories: emotional abuse (e.g., 'family members called me stupid, lazy, or ugly'), physical abuse (e.g., 'someone from my family hit me so hard that I had to see a doctor or go to the hospital'), sexual abuse (e.g., 'someone tried to touch me in a sexual way or made me touch him/her in a sexual way'), emotional neglect ('e.g., I felt loved'), physical neglect (e.g., 'someone took me to the doctor when it was necessary').Each category contains five items.Items 2,5,7,13,19,26, and 28 were reversed before taking the sum score for current analyses.If data was missing for one time point but was available for other timepoints, the missing datapoint was imputed for that subject based on the average of the remaining time points.Imputation was done within the five categories separately (sessions / % per category: emotional abuse: 15/0.0054,physical abuse: 4/0.0014, sexual abuse: 4/0.0014, emotional neglect: 17/0.0061,physical neglect: 19/0.0068).Moreover, the CTQ was part of the questionnaires given only to participants of the U-change cohort, meaning that the measurement time point did not match that of other questionnaires.As CTQ items are also less timepoint specific compared to e.g., the LEQ, we used the average CTQ rating across sessions for our analyses.
The Alabama Parenting Questionnaire (APQ) as included in the NSPN study contains 15 items asking about parenting styles.These 15 items are a combination of 9 items from the original version of the APQ 44 , the Corporal Punishment scale (3 items), and the Involvement scale (3 items).Participants rated the frequency of occurrences of certain parenting styles in their family on a five-point scale ('never' to 'always), asking about positive parenting (3 items), inconsistent discipline (3), poor supervision (3 items), involvement (3 items), and corporal punishment (3 items).We reversed the 'positive parenting' and 'involvement' items so that higher scores reflect more adverse parenting styles.Imputation was performed in the same way as described for the CTQ.Imputed sessions/% per category: positive parenting: 30/0.0054,inconsistent parenting: 51/0.0062,poor supervision: 48/0.0086,involvement: 50/0.0090,corporal punishment: 12/0.0022.
The Measure of Parenting Style (MOPS; 45 ) measures dysfunctional parenting practices in 30 items, 15 for each parent.The questionnaire covers indifference/neglect (6 items, e.g., 'Was uninterested in me'), over-control (4 items, e.g., 'Sought to make me feel guilty'), and abuse (5 items, e.g.'Physically violent or abusive of me ').Higher scores reflect more dysfunctional parenting styles in all sub-scales.Scores were first summed within categories and then across categories for the current analyses.
Socioeconomic status was indirectly derived from the Index of Multiple Deprivation (IDM) associated with the area a participant lived in.The MDI is based on regional income, employment rate, education, health and health service, crime rates, barriers to housing and serviced, living environment and others.IDMs were imputed based on the mean if missing (15 sessions = 0.01%).

Cross-sectional effects of ResPSF and ∆ ResPSF on MT
In order to assess whether the observed positive association between ∆ResPSF and ∆MT in the prefrontal cortex co-occurs with 1) pre-existing hypo-or hypermyelination, or 2) cross-sectional differences in baseline MT as a correlate of ResPSF, we applied two general linear models to baseline MT data:

Supplementary Figure S2 .
Representativeness of the MRI subsample.The MRI subsample comprises the n = 141 individuals included in the main analyses linking changes in resilient psychosocial functioning (ResPSF) scores to myeloarchitectonic and functional maturation.The Non-MRI sample comprises all other individuals that were included in behavioral analyses only and for whom respective behavioral data were available: In A), the non-MRI sub-sample includes n = 314 individuals with longitudinal ∆ResPSF scores and n=885 individuals with distress scores.In B), the non-MRI sub-sample comprises n = 1457 individuals with all risk exposure assessments completed.Distributions show density plots.APQ = Alabama Parenting Questionnaire; MOPS = Measure of Parenting Style; LEQ = Life Events Questionnaire; IMD = Index of Mean Deprivation.

Table S1 .
Factor loadings of questionnaires included for the computation of a general distress factor.

Table S2 .
Distribution of distress and resilient psychosocial functioning (ResPSF) scores with respect to sex (mean +/-SD) and age in the prediction sample (i.e., 712 individuals included for the computation of ResPSF; 455 with repeated measures) and the imaging sample (i.e., 141 individuals for which both imaging and behavioral data was available for two time points).